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We present a theory of electron-mediated interaction between adatoms in graphene. In the case of 
resonant scattering, relevant for hydrogentated graphene, a long-range 1/r interaction is found. This 
interaction can be viewed as a fermionic analog of the Casimir interaction, in which massless fermions 
play the role of photons. The interaction is an attraction or a repulsion depending on whether the 
adatoms reside on the same sublattice or on different sublattices, with attraction dominating for 
adatoms randomly distributed over both sublattices. The attractive nature of these forces creates 
an instability under which adatoms tend to aggregate. 

PACS numbers: 



Unique transport characteristics of graphene make it a 
strong candidate for replacing silicon in future electronic 
devices [l| . Functionalizing graphene by controllable ox- 
idation 0, S| or hydrogenation [j, 0] can change its prop- 
erties in new, unexpected ways. In particular, when hy- 
drogen adatoms bind to graphene, the orbital state of 
each functionalized carbon atom changes from sp^ to sp^ 
configuration, removing 7r-electrons from the conduction 
band and turning graphene into a semiconductor . Re- 
markably, the metallic properties can be fully recovered 
after dehydrogenation This provides a unique tool to 
control electronic properties of this material [7|, iS] • 

One of interesting questions posed by the experiment 
has to do with the interaction between adatoms me- 
diated by electron scattering. As we shall see, resonant 
scattering on the midgap states localized on adatoms 
0, 0, El leads to dramatic enhancement of interac- 
tion, making it long-ranged. We find that the interac- 
tion energy falls off very slowly, approximately as inverse 
distance between the adatoms, U{r) ^ r^^ . The sign of 
interaction depends on the sublattice type: two atoms re- 
siding on different sublattices {A and B) attract, whereas 
atoms on the same sublattice repel [see Eqs. ([TT|) . (fT3|) ]. 

The interaction is stronger than the long-range in- 



teraction between adatoms on surfaces of metals [12,112, 
lJ|, which is of a Friedel-oscillation (FO) character. The 
FO interaction falls off as when it is mediated by elec- 
tronic states on the surface, and as , when mediated 
by the states in the bulk 15, [iBl- The FO interaction 



can occur in graphene [l3|- Long-range interaction can 
lead to fascinating collective behavior of adatoms, such 
as self-organization into chains 14 1 and superlattices [3| ■ 
The interaction analyzed in this work can be inter- 
preted as a fermionic Casimir effect. The Casimir inter- 
action between two bodies (or, atoms) arises due to scat- 
tering of virtual photons. For each of the bodies, in the 
presence of the second body, angular distribution of the 
flux of incident virtual photons is somewhat anisotropic, 
giving rise to a net attraction force. This interaction is of 
a generic character (fermionic Casimir effect was recently 
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FIG. 1: Electron-mediated interaction between adatoms in 
graphene modeled by a hard-core potential: numerical results 
(black dots) and mean field theory, Ea. (|18p (red line). The net 
interaction is a repulsion when adatoms are randomly placed 
on one of the sublattices, and an attraction when they are 
equally distributed over both sublattices. The 3/2 power law 
U]) provides an accurate fit to the numerical results with the 
best fit values eo = — 0.75eV (top curve) and £o = 1.3 eV 
(bottom curve). System of size 48 x 82 was used for simu- 
lation, each data point was averaged over 20 realizations of 
randomly generated adatom configurations. Inset: Attracting 
adatoms tend to aggregate. Phase diagram obtained from the 
free energy (fT9)l is shown. 



analyzed in one-dimensional systems [l9|). 

We find that attraction between atoms on different 
sublattices is stronger by a logarithmic factor than re- 
pulsion within the same sublattice. The net interaction 
of atoms equally distributed among the two sublattices 
is thus an attraction, characterized by the energy density 



E{n) 



-Son 



3/2 



1.3 eV, 



(1) 



per carbon atom (see FiglT]), where n is the fraction of 
hydrogenated carbon atoms. The prefactor in ([T]) may 
have a weak logarithmic dependence on n. 

We emphasize that the interaction energy in this case 
cannot be treated as a sum of pairwise two-particle inter- 
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FIG. 2: Electron-mediated interaction between adatoms de- 
pends on the type of sublattice: atoms on different sub- 
lattices, A and B, attract (a), whereas atoms on the same 
sublattice repel (b) [see Eg. pip and Ea. (|13p ]. The interac- 
tion is modulated by a prefactor which takes different val- 
ues on the three sub-sublattices marked by 2, 2' and 2"; 
(a) |sin(Kr4-(/.)| = |sin(^|, |sin(0-F ^)|, | sin(</. - ^)|; (b) 
cos^(Kr) = 1, J, i. The modulation results from interference 
between electronic states in valleys K and K' . 



actions (indeed, summing 1/r interactions over the entire 
space would give a divergence). The situation resembles 
that of Casimir forces, which are of an essentially non- 
pairwise nature. To treat the interaction mediated by 
electrons one must account for the change in electronic 
states at the energies e < hvon^^"^, resulting from elec- 
tron scattering on the adatoms {vq is the electron Fermi 
velocity) . This leads to interaction energy per adatom of 
order hv^/r, with r = •nT^I'^ the typical distance between 
adatoms, in agreement with r?!"^ scaling, Eq.(IT]). 

Attraction can lead to instability of a homogeneous 
phase and adatom aggregation. Characteristic time 
scales for such processes are controled by the rates of ad- 
sorption and desorption, or diffusion, whichever is faster. 
Compression of the graphene lattice, resulting from at- 
traction between adatoms, may expalain the observed re- 
duction of the lattice constant [1] , which is at odds with 
first-principles calculations . 

Interaction between hydrogen adatoms could also re- 
sult from corrugation of the graphene sheet caused by the 
stress around tetrahedral s-p^ bonds. Numerical evidence 
suggests, however, that such corrugation is limited to the 
range of at most a few lattice constants [8], rendering this 
type of interaction effectively short-ranged. 

The problem of electrons scattering on impurities can 
be described by a tight-binding Hamiltonian 

if = ^(tk?AL^V'k.S+h.C.)+ M„(x)V'i,„^x,a. (2) 

Here u^(b)(x) is adatoms' potential on sublattices A(B), 
andtk = io(l+e-**'^i+e-*'=^), with to ~ 3.1 eV the hop- 
ping amplitude and ei(2) the basis vectors (see Figllb). 

The interaction between adatoms can be conveniently 
analyzed in the Matsubara Greens function framework 
using the thermodynamical potential = Tr InG 



[20|. For two adatoms, we write — Gq — Vi(x — 
Xi ) — V2 (x— X2 ) . Resumming perturbation series in terms 
of the T-matrices of each adatom, we obtain 



^ = Tr In (1 - T1G12T2G21] 



(3) 



Here 6*12 is the free-particle Greens function in position 
representation, evaluated between the points Xi and X2 
(similar representation was used recently in a study of 
Casimir forces 21, 2^). 

In this paper we shall use the s-wave resonant scatter- 
ing approximation, 



Toiie) 



ieh-i{W/\e\) + 6' 



\d\ « ly « 3to, (4) 



as appropriate for short-range scatterers at low ener- 
gies. Here W is the electron half-bandwidth, and the 
parameter S describes detuning of resonance from the 
Dirac point. An expression of this form can be ob- 
tained for a delta-function potential, m(x) = V5{'k — 
xi), in which case the T-matrix is given by T{e) = 

V/ (1 + -^ie In ^) [i, K El- For hydrogen adsorbed 
on graphene, the presence of a resonance peak close to 
the Dirac point, Eq.(l4|), was confirmed by first-principles 
calculations Q. 

The real-space Greens function can be written as 



G{ie,r) = - 



zkr 



(27r)2 £2 



ie tk 

tl ie 



(5) 



where the matrix accounts for the A and B sublattices. 
The Greens function takes on different form for the end 
points on different sublattices: 



G(ie,r) 



Gaa Gab 
Gba Gbb 



(6) 



In the low-energy approximation we expand in the 
vicinity of points K, K' = — K to obtain « vq{'^Px — 
ipy), where p = k =p K is the momentum relative to the 
K (K') point, and Vq — ^Iq is the Fermi velocity. Adding 
contributions of the states near K and K', we obtain 

^ ^ iecos(Kr) , r 

Gaa ^ Gbb ^ '-^Ko{er), r = — , (7) 



vo 



Gab = 



£ sin(Kr 



-Ki (ef) 



(8) 



where cj) is the angle between r and K (see Figl2^), and 
Kqa denote modified Bessel functions of the second kind, 

r(i/+i)2'' foo cos ztdt 



K,{z) 



Jo 



The function Gba can 



W (l + t2)''+l/2 ■ 

be obatined from the relation Gba{j^) = G'^^(— r), giving 



£ sin(Kr - ^) 

Gba = 2 -^1 i'^^) ■ 



(9) 
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We first consider two adatoms on different sublattices 
(see Figl2^). At distances r ^ Hvq/T, approximating the 
sum in = -TY.e„ 1^ (l - T^(te)G\j,{ie^,v)) by an 
integral / and using Eq.®, we find 



^12 



de 
2^ 



log 1 



£2 siT?{K.Y + (j))Kl{ef) 
{ie\og{W/e)+5f 



(10) 



This result further simplifies for relatively short distances 
r ^ Hvq/S. The integral can be evaluated using the 
asymptotic formula Ki{x <C 1) ~ 1/x and replacing 
hi{W/e) by ln{rW/hvQ) with logarithmic accuracy. Set- 
ting 5 = and using the identity dx\n{l + u/x^) ~ 
-K^Ju we integrate over e to obtain 



ftwol sin(Kr + 0)| 
r log(r/a) 



(11) 

where a = Hvq/W. The interaction has a negative sign, 
corresponding to attraction of adatoms. 

Interestingly, due to the factor | sin(Kr + cj))\ in the 
above expression, the interaction oscillates on the lattice 
scale. This oscillation results from interference of the 
contributions due to fermions from K and K' valleys. 

The meaning of the factor | sin(Kr + ip)\ can be seen 
more clearly by considering it separately on each of the 
three sub-sublattices, which have period \/3 times the 
period of the A or B sublattice (see Figl2^). Since e*^"" 
takes values 1, e^in/a^ ^j^fj g^'^*/^, the same on each of the 
three sub-sublattices, the angular dependence in Eq. (fTT|) 
is given by | sin((/))|, | sm{(j) + 2n/3)\, or | sin(0 -I- 47r/3)| 
in each of the three cases. 

For a pair of adatoms residing on the same sub- 
lattice {A or B), the interaction energy is ^li2 — 
-TEs^ ln(l- (To(i£„)GA^(ie„,r))2), giving 

£2 cos^(Kr)K^{ef) 



de 
2^ 



In 1 



{ie\n{W/e) +Sy 



(12) 



We note a different sign under the log in this expres- 
sion as compared to Eq. pil| . which arises because Gaa 
is imaginary-valued, whereas Gab is real-valued. The in- 
tegral over £ is dominated by the region (5 ^ |£| ^ vq/t, 
since -^^0(2;) decreases exponentially for x > 1. For such 
£, and for ln(T/Fr/uo) > 1, the ratio KQ{ef)/\n{W/s) is 
small in most of the integration domain [Kq{x 1) ~ 
— logo;]. Thus we can Taylor-expand the log and, with 
logarithmic accuracy, integrate over e using the identity 
Kl{x)dx — 7r2/4, to obtain 



4rlog (r/a 



— cos2(Kr). (13) 



UAAia < r < hvo/S) 



The factor cos2(Kr) in Eq.([T3]), describing interference 
between two valleys, takes constant value on each of the 
three sub-sublattices with period \/3 (see Figl^jD). Ana- 
lyzing it as above we find that cos^ (Kr) = 1 for adatoms 
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FIG. 3: The interaction (|13p . is sensitive to the resonant 
form of the T- matrix, Eq.Q. For a nonzero detuning 5, the 
interaction retains the 1/r form at distances r < hvo/&, de- 
creasing more rapidly at larger r. When the system is doped 
away from neutrality, similar behavior is expected at distances 
shorter than the Fermi wavelength, r < Ap. 



residing on the same sub-sublattice, and cos2(Kr) = 1/4 
when adatoms reside on different sub-sublattices. 

The energy of interaction for adatoms on the same sub- 
lattice, Eg. lfT^ . is positive, which means that in this case 
adatoms repel each other. This repulsion is logarithmi- 
cally weaker than the attraction found for atoms on dif- 
ferent sublattices, Eg. pT]) . We thus expect the net in- 
teraction for a system of many adatoms randomly placed 
on both sublattices to be dominated by attraction. 

The repulsion (fT3)) will be greatest for the next-nearest 
carbon atoms. Interestingly, in an STM experiment [2^ 
it was found that chemisorbed hydrogen atoms can reside 
on the nearest or next-next-nearest sites of the carbon 
lattice, but never on the next-nearest sites. This behavior 
is consistent with our results, Eqs. p3)) . pT|l . 

Next, we analyze interaction in a system of adatoms 
at a finite concentration. Since electronic states with 
wavelengths exceeding the distance between adatoms, 
A ^ d = n~^l'^ ^ are strongly perturbed by scattering, 
this interaction is of non-pairwise character. For rela- 
tively high densities, n > 6'^ /fi^VQ, the interaction can be 
estimated using the results for 6 = 0. This gives an en- 
ergy £* = hvo/d per adatom, leading to the n^/^ scaling 
for the energy density vs. adatom concentration, Eq.(IT]). 

This behavior was confirmed by numerical analysis of 
the tight-binding problem whereby adatoms were 
modeled by a local potential taking values exceeding tg. 
Given a random configuration of N adatoms, we diago- 
nalize the Hamiltonian and sum all negative eigenvalues 
to evaluate the total energy, E{N) = J2e <o^"- '^^^ 
dependence on A'^ is dominated by a contribution lin- 
ear in N, E{N) = Eo + AqN + AiN^/^, which repre- 
sents a chemical potential of an adatom. Subtracting the 
linear part, we recover the interaction AE{N) cx N^^'^ 
(see Fig.[T]). Alternatively, one can choose to evaluate 
E{N) as a sum over the lower half of the spectrum. This 
changes somewhat the linear term, leaving the N^^'^ con- 
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tribution essentially the same. 

The sign of interaction is that of attraction when 
adatoms are evenly spread over both sublattices. In this 
case, the best-fit value of the prefactor in the scaling re- 
lation HI) is found to be eo « 0.42 to- With to = 3.1 eV 
this gives eo ~ 1-3 eV. In contrast, when all adatoms are 
placed on one sublattice, a repulsive interaction is found, 
eo « —0.24^0 = — 0.75eV. This is in agreement with the 
signs of pairwise interaction discussed above. 

To test these numerical results against an analytic ap- 
proach, we use disorder-averaged Greens functions found 
in a selfconsistent mean-field approximation, in which 
point-like adatoms are replaced by a constant field: 



G ^{ie, k) — is — 
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(14) 



where ni = 2nl??l'^c? is adatoms' density per sublattice, 
a = 1.42 A is carbon spacing. Solving the selfconsistency 
condition ([M]) with logarithmic accuracy, we find 



sen e 
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^=2 ■ —V 4 
The energy density of the system can be written as 



A^ln-^ = 7ruo"-i- (15) 
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27ri' 
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z - en 



de 
2^ 



ieTrG(ze) (16) 



giving the critical temperature = eo/2\/3 « 4200 K. 
Since temperature during hydrogenation is substantially 
below [5|, the adatoms are expected to self-organize 
into high and low-density droplets. 

Even if spatial diffusion of hydrogen is slow, as may be 
the case in initial stages of self-organization termi- 
nated by freezing in a low-temperature state would result 
in macroscopic inhomogeneities. Such inhomogeneities of 
the hydrogenated state were indeed observed in the TEM 
diffraction images described in Ref. [3| • It was also noted 
that dehydrogenation restores homogeneity, pointing to 
an intrinsic character of this effect. 

The attraction between "frozen" adatoms would cre- 
ate a lateral stress. Treating the occupancy n as strain- 
independent, we have 



a = -dE{n)/d\nV ^ ^ko|?^^/^ 



(20) 



where an empirical relation dto/da « —to /a is used to 
describe the change in tQ. Such stress would lead to com- 
pression of the graphene matrix. This is consistent with 
the reduction in lattice period upon hydrogenation ob- 
served in experiment ^Ss] . 

We are grateful to A. K. Gcim, R. L. Jaffe, and 
K. S. Novoselov for useful discussions. 



where Sa is the spectrum, and the contour integral is 
taken over the imaginary axis and a half-circle at infin- 
ity. The trace of G is identical to that in the self-energy of 
a T-matrix, giving TrG'(ie) = — 2ieln(M^/|e|)/7rwo. Sub- 
tracting the contribution due to free Dirac fermions, we 
obtain the change in total energy due to adatoms. 



de 
(ttw)^ 



e In 



W 



e In 
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(17) 



The function under the integral is even, positive, and 
approximately constant at |e| A, taking on a value 
proportional to n (with logarithmic corrections). At 
< e ^ A the function is increasing from zero to the 
asymptotic value at large e. This behavior is in agree- 
ment with expectation of a leading contribution 5E cx n 
and a negative r?^"^ part describing interaction. Subtract- 
ing the part linear in n, and dividing by the density of 
carbon atoms 7io, we find the interaction energy 



AK, 



8A3 



In 
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no 



33/2^2 ' 



(18) 

per carbon atom. This formula agrees very well with our 
numerical results (see red curve in Fig[T]). 

A long-range attraction between adatoms can drive 
thermodynamic instability. This can be seen most eas- 
ily from phase diagram, obtained from the free energy 
F = E{n) — TS{n) (see FigH] inset). In our case, 

F = -eon^''^ + T {n\mi+ [l - n)\n{l - n)) , (19) 
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